
    
    DREst = function(data, cnt,alpha.pmle,beta.pmle,gamma,y.hat,l1.hat, SEED = 1, dr.warm=NULL,dr.warm.obj=NULL){
        # Output: ("est","nll","time", "contrast")
        # nll is just filled with 0---DR has nothing to do with nll
        
        ptm = proc.time()
        point.est <- dr.estimate.noiterate(data, cnt, alpha.pmle,beta.pmle, gamma,y.hat=y.hat
                                           ,l1.hat=l1.hat, Hessian = F, SEED=SEED, dr.warm=dr.warm,dr.warm.obj=dr.warm.obj)
        time = (proc.time() - ptm)[3] 
        
        
        alpha     = point.est
        beta      = beta.pmle
        est   = abind(alpha,beta, along=3)
        
        time = array(time, dim(est))
        # GetContrast(data, alpha.dr, beta.dr) #DEBUG
        contrast = GetContrast(data, alpha, beta)
        contrast = array(contrast, dim(est))
        # contrast = array(0, dim(est))
        
        nll  = array(0, dim(est))
        output = abind(est, nll, time, contrast, along = 4)

        return(output)
        
    }
    
    DREst.Iter = function(data, cnt,alpha.pmle,beta.pmle,gamma, max.step = 500, thres = 1e-4){
        # Output: ("est","nll","time", "contrast")
        # nll is just filled with 0---DR has nothing to do with nll
        Optim = function(par.update, fn, parn, pars, data,cnt){
            optim(par.update, fn, parn=parn, pars = pars, data = data,cnt=cnt, 
                  control=list(maxit=max.step), method =  "L-BFGS-B",
                  lower = -rep(10,10), upper = rep(10,10))
        } 
        
        pars = alpha.pmle
        beta = beta.pmle
        gamma = gamma
        
        ## Optimization 
        ptm = proc.time()
        Diff = function(x,y) sum((x-y)^2)/sum(x^2+thres)
        diff = thres + 1; step = 0
        while(diff > thres & step < max.step){
            if(step %% 10 == 0) cat("This is the ", step, "th step. Diff = ", diff, "\n")
            step = step + 1
            opt1 = Optim(pars[1,], dr.objective.iter, "alpha1", pars, data,cnt)
            diff1 = Diff(opt1$par,pars[1,])
            pars[1,] = opt1$par
            
            opt2 = Optim(pars[2:5,],dr.objective.iter, "alpha2345", pars, data,cnt)
            diff2 = Diff(opt2$par,pars[2:5,])
            pars[2:5,] = opt2$par
            
            diff = max(diff1, diff2)
            
        }
        time = (proc.time() - ptm)[3] 
        
        
        alpha     = pars
        beta      = beta.pmle
        est   = abind(alpha,beta, along=3)
        
        time = array(time, dim(est))
        
        contrast = GetContrast(data, alpha, beta)
        contrast = array(contrast, dim(est))
        
        nll  = array(0, dim(est))
        output = abind(est, nll, time, contrast, along = 4)
        
        return(output)
        
    }
    
 